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Structural, vibrational and thermal properties of densified sodium silicate (NS2) are investigated 
with classical molecular dynamics simulations of the glass and the liquid state. A systematic investi- 
gation of the glass structure with respect to density was performed. We observe a repolymerization 
of the network manifested by a transition from a tetrahedral to an octahedral silicon environment, 
the decrease of the amount of non-bridging oxygen atoms and the appearance of three-fold coor- 
dinated oxygen atoms (triclusters). Anomalous changes in the medium range order are observed, 
the first sharp diffraction peak showing a minimum of its full-width at half maximum according to 
density. The previously reported vibrational trends in densified glasses are observed, such as the 
shift of the Boson peak intensity to higher frequencies and the decrease of its intensity. Finally, we 
show that the thermal behavior of the liquid can be reproduced by the Birch-Murnaghan equation 
of states, thus allowing us to compute the isothermal compressibility. 
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I. INTRODUCTION 

In the field of oxides, silicate glasses and melts have re- 
ceived a huge attention for their important applications 
in materials science and geophysics, such as magma dy- 
namics and properties. Pressure (or density) is obviously 
one of the most important thermodynamic variable for 
geochemical processes in the mantle and crust. Indeed, 
interesting macroscopic properties of silicate melts, such 
as viscosity or diffusion, show significant changes with 
pressure*^ 

Many experimental studies on silicate glasses, the 
base material for various multi-components silicate sys- 
tems, have suggested that those macroscopic proper- 
ties were related to atomic-scale structural changes 3 
such as angles^ or coordination number— ~—. Densified 
sodium silicate is a very interesting system to be in- 
vestigated as it shows the effect of polymerization and 
depolymerization 5 ' 7 . Indeed, in the silica network, Si 
tetrahedrons are connected by bridging oxygen atoms 
(BOs). Sodium silicate is usually described as a base sil- 
ica network which is depolymerized by the sodium atoms. 
In this view, sodium cations break Si-BO-Si bonds and 
induce non-bringding oxygen atoms (NBOs). On the con- 
trary, pressure tends to repolymerize the network by a 
global increase of coordination numbers. 

Sodium silicate glass has already been extensively 
studied at ambient pressure using Molecular Dynamics 
(MD). The first reported MD simulation of sodium sili- 
cate glass in 1979 was based on a very small system (200 
atoms) but it is remarkable to see that it presented a 
very reasonable structural description of the glass. Since 
this work, the used potentials have been continuously 
improved to get a better reproduction of experimental 
results. Progress in computing facilities progressively al- 
lowed to reach longer time scales, thus making it possible 
to study diffusion at lower temperature and rheological 
properties^. The possibility to simulate larger systems 
has also permitted to put in evidence the existence of 



inhomogeneities and preferential diffusion pathways for 
sodium atomsiSr— . Simulations from Cormack and co- 
workers^— have shown a very good agreement with ex- 
perimental results on structure. Vibrational 11 ' 20 ! 21 and 
elastic 22 properties of the glass at ambient pressure have 
also been studied and successfully compared to experi- 
mental data. Using superomputers, large scale classical 
simulations have recently been performed 23 , as well as ab 
initio Molecular Dynamics simulations 2 ^—. However, to 
our knowledge, no systematic study of the evolution of 
the system according to pressure has been performed so 
far. 

We present here Molecular Dynamics simulation allow- 
ing a systematic description of the structural, vibrational 
and thermodynamics properties of densified glassy and 
liquid sodium silicate. We focus in one particular com- 
position (NS2) and study the properties with increasing 
density. Results show that a transition from tetrahedral 
to octahedral silicon environment occurs and that the 
medium range order shows anomalous changes. Vibra- 
tional properties are also found to be very sensitive to 
pressure and we report some trends about the behavior 
of the Boson peak according to density. Eventually, an 
equation of state model is proposed, thus allowing the 
computation of the isothermal compressibility. 

The article is organized as follows. In section II, we 
present the numerical model and methodology that has 
been used. In section III, we report structural, topo- 
logical and vibrational results of the glass. In section 
IV, thermodynamics and structural results of the liquid 
state are presented. Finally, section V summarizes these 
results. 



II. SIMULATION DETAILS 

As just mentioned, (Na20) K - (Si02)i-a; with x=0.30 
system has been chosen (close to the so-called NS2 system 
with x=0.33). The simulated system is composed of N = 



3000 atoms (700 Si, 1700 O and 600 Na), placed in a cubic 
box of various lengths L to study different densities (from 
1.5 g/cm 3 to 5.4 g/cm 3 ). To do so, all the simulations 
were run in the canonical ensemble (NVT). The room 
temperature density^ 7 - of 2.466 g/cm 3 is obtained with 
L=34.43A. A computed pressure P = -1.6 GPa is found 
in the glass at this density. 

To take into account the oxidation state of atoms^ 8 -, 
partial charges are used for the Coulomb interaction, 
while the short-range Buckingham potential is of the 
form : 



V l3 {r) = A l3 exp{~— ) - ^ 
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where Aij , gy and Cy are parameters which have been 
fitted by Teter 18 . Usually, the Buckingham potential can 
induce spurious effects at high temperature (as V(r) can 
go to negative infinity when r is close to zero, which 
leads to a collapse of the interacting atoms 28 . As de- 
scribed ir>i£, a repulsive term Bij /r nij was introduced at 
short distance in order for the potential energy and its 
derivative to be continuous at r$ to avoid this issue. 

This potential has been extensively used by Cormack 
et a L 17 ' 18 and has revealed a very good description of the 
glass at room density for various compositions. The ef- 
fect of pressure on such systems using classical Molecular 
Dynamics has been considered 2 — only at high tempera- 
ture for the NS4 silicate by using a Born-Mayer interac- 
tion potential fairly similar to the one that is presently 
used. While, at ambient pressure, the use of a Coulomb 
interaction with fixed partial charges is supported by the 
ionic character of the interactions and the absence of 
charge transfer, one may wonder to what extent fixed 
charges can be still considered with increasing pressure. 
While we are not aware of any report of densified sili- 
cates, a recent ab initio Molecular Dynamics study (in 
which electrons and charge transfer are explicitly com- 
puted) on an oxide network-forming glass under high 
pressur e) 30 ' 31 has not shown any deformations of the elec- 
tronic cloud that would be significant enough for ambient 
pressure pseudopotentials to be modified. The mentioned 
exampl e 30 ' 31 , the consistency of the presently reported 
results and the fact that the present potential was suc- 
cessfully used to reproduce a diffusion anomaly of O and 
Si atoms with increasing densit y 16 ' 32 , also observed in 
pure silica 33 or water 34 , suggest that a certain degree of 
confidence can be expected. 

Classical Molecular Dynamics simulations were per- 
formed using the DLPOLY package 3 -. The equations 
of motion were integrated with the Verlet-Leapfrog algo- 
rithm, using a timestep of 2.0 fs. Coulomb interactions 
were evaluated by the Ewald summation method with a 
cutoff of 12.0 A. The short-range interaction cutoff was 
chosen at 8.0 A. As mentioned, the simulations were 
run in the canonical ensemble (NVT) with a Bcrendsen 
thermostat. 

For each density, the system was first equilibrated at 
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FIG. 1: (Color online) Total radial correlation function of MD 
modeled sodium silicate glasses for increasing densities and 
comparison with neutron diffraction studies (white rounds) 
from the work of Wright et al.— (Neutron diffraction data). 



6000 K during 10 6 steps (2ns). Each melt was then con- 
tinuously cooled down to the selected temperature (from 
300 K to 4000 K) using a cooling rate of 10 K/ps. 



III. GLASS 
A. Real space properties 

1. Total radial correlation functions 

The total correlation functions gj(r) for increasing 
densities are shown in Fig. Q] To check the validity of 
the simulated glass, comparison with experimental data 
(neutron diffraction from the work of Wright et al. 36 ) 
at room pressure was made. We recover the same level 
of agreement than in previous studie d 17 ' 19 . However, we 
notice an increased structured system with main peaks 
being sharper as compared to experiments. This com- 
parison has also been done by Cormacki 7 -. Using the 
same potential, a better agreement has been observed 
by broadening the total correlation functions^ 7 - The po- 
sition of the first Si-0 peak is well reproduced, but is 



found to be sharper than in experiments. The position 
of the second 0-0 peak is also well reproduced, suggest- 
ing a realistic O-Si-0 angle in simulation. On the other 
hand, simulation produces a peak at 3.1 A arising from 
Si-Si correlations (see below) which is not present in ex- 
periments but merged with other contributions in the 
region 3-4A. It means that the inter-tetrahedral angle 
Si-O-Si may be underestimated with respect to experi- 
ments. This angle turns out to be highly sensitive to 
the employed potential. A detailed discussion about the 
ability of the different potentials to reproduce the Si-O-Si 
angle can be found mi£. 

As density increases, the first Si-0 peak does not show 
any shift in position but becomes broader, suggesting 
an increased disorder in the network, manifested by in- 
creased coordination numbers (integral of the first peak) . 
As observed on the partial gi{r) distributions (see below), 
the second peak is shifted to lower r and becomes broader. 



2. Partial radial correlation functions 



The partial radial correlation functions g»(r) have been 
computed from the pair correlation functions gy (r) : 



9i{r) 



1 n 
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(2) 
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We have split the analysis according to BO and NBO. 
These functions are shown in Fig. [2] for increasing den- 
sities. While the position of the first peak in gs\ (Si-0 
correlations) does not show any significant change, an in- 
crease in the shoulder on the lower r side of the second 
peak (Si-Si correlations) is observed as density increases, 
suggesting that the Si-O-Si angle decreases. As men- 
tioned previously, the environment of the BO and NBO 
are studied separately using c/bo and <7nbo- For both, 
the position of the first peak (O-Si correlations) remains 
the same, but important changes take place with density 
change for the second-neighbor correlation. The second 
peak (0-0 correlations) is shifted to lower r and the dis- 
tribution becomes broader. In the (7nbo partial correla- 
tion function, one notices the growth of a peak (at 2.4A 
for g = 3.5 g/cm 3 ) which contributes only to a shoulder 
of the main peak at 2.6A for g = 2.5 g/cm 3 . This also 
suggests that densification affects the O-Si-0 and Si-O- 
Si angles rather than the Si-0 distance between nearest 
neighbors. Finally, the first peak of gNaOOj associated 
with Na-0 correlations, is shifted to lower r. The de- 
crease in the Na-0 distance with pressure has also been 
observed using NMR by Lee£. The Na-centered pair dis- 
tribution functions are highly sensitive to density change 
and this is not surprising as it involves non-directional 
bonds. However, we notice that the increase of density 
leads to a better defined first peak whose height increases 
with the density. 



3. Coordination numbers 

In pure silica, the network in fully connected and the 
coordination number CN of Si and O atoms are found 
to be 4 and 2, in agreement with the stoichiometry of 
the glass (CN Si N S i = CN N ). This is not the case 
in sodium silicates since Na atoms create NBOs, thus 
disrupting the network. 

The distributions of IV, V and Vl-fold coordinated sil- 
icon atoms (Si IV , Si v and Si VI ) can be obtained by enu- 
merating the number of oxygen neighbors in the first co- 
ordination shell of each silicon atom. These populations 
are shown in Fig. [3^ for each CN. The fraction of tetra- 
hedral Si IV atoms starts to drop from g = 2.7 g/cm 3 (P 
~ 1 GPa). At the same density, the fraction of Si v atoms 
grows and reaches a maximum around g = 4.0 g/cm 3 (P 
~ 28 GPa) prior to a continuous decrease as density in- 
creases. The fraction of octahedral Si VI atoms increases 
from g — 3.1 g/cm 3 (P = 5 GPa) and this basic structure 
becomes predominant at high density. These trends are 
rather usual in densified silicates. In amorphous silica, 
simulations from Tse 3 ^ predicted the increase of the Si 
CN to 5 at 15 GPa and up to 6 at 20 GPa. That trend 
was confirmed by simulations from Horbach 39 . The ap- 
pearance of Si v and Si VI in densified sodium silicate has 
been confirmed experimentally using NMRj^ 

The environment of oxygen atoms has been analyzed 
in the same fashion, i.e. by enumerating the number of 
silicon atoms in the first coordination shell of each oxy- 
gen atom. Here, Na atoms are not taken into account, 
this in order to distinguish BO from NBO and thus to 
split the Si CN analysis from the one involving the Q™ 
speciation which will be detailed below. Thus, O 1 refers 
to the oxygen atoms that are surrounded by only one sili- 
con atom (i.e. NBO atoms). At low density, the fraction 
of NBO can be determined by x, the amount of soda, 
as each sodium atom creates one NBO. The fraction of 
NBO f N BO is thus given by f N BO = N NB o/No = 2x/(2- 
x). At x = 0.3, fNBO ~ 0.35, which is consistent with 
simulation results in Fig. (3Jd. The fraction of O 1 drops 
for densities larger than g = 2.6 g/cm 3 . At this density, 
the fraction of O 11 starts to increase, reaches a maximum 
at g — 3.8 g/cm 3 and decreases at higher density. The 
present findings clearly indicate a repolymerization of the 
network through the creation of density induced Si-BO- 
Si connections at the expense of Si-NBO ones. They 
are consistent with the decrease of the fraction of NBO 
found experimentally from NMR in densified silicates^. 
As a consequence, the model of the network modifier Na 
atom simply given by stoichiometry (one Na atom involv- 
ing the appearance of one NBO atom) does not remain 
valid for g > 3 g/cm 3 . Ultimately, O m are found at high 
density and their fraction grows up to nearly 50% at g 
= 5.5 g/cm 3 . Note that 3-fold O atoms (termed triclus- 
ters) have already been found both in experiments and 
in simulations, for example in aluminosilicate glasses^ . 

At g = 5.5 g/cm 3 , the fraction of NBOs is very low 
(~ 3%) so that the Si/O network can be considered as 
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FIG. 2: (Color online) Partial radial correlation function gsi(r) 
2.5, 3.5, 4.5 g/cm3. 



<?Na(r), 9Bo{r) and <?nbo(?") at different selected densities g = 



being fully connected as in pure silica, thus allowing us 
to check the agreement between the stoichiometry of the 
system (Si02.43) and the computed coordination num- 
bers. On average, we find CNsi = 5.90 and CNo = 2.43 
so that the stoichiometry of the glass is satisfied (CNsiNsi 
« CNqNo). 



These results show that the network undergoes strong 
topological changes as density increases. The initial 
tetrahedral silicon environment becomes octahedral at 
high density, consistently with the decrease of the O-Si-0 
angle (see below). On the other hand, a transition from 
2-fold to 3-fold coordinated oxygen atoms is observed, 
which is once again consistent with the decrease of the 
Si-BO-Si angle (see below). 



4- Q n populations 

As mentioned earlier, changes in the glass network can 
also be characterized by the Q" distribution analysis. We 
remind that a Q n silicon atom is defined as an atom 
linked with n bridging oxygen atoms. Defining BO and 
NBO at high density needs a careful analysis since 3-fold 
coordinated oxygen atoms can be found. NBOs are here 
thus defined as oxygen atoms connected to only one Si. 
BOs are defined as oxygen atoms that are not NBOs. 

At ambient pressure, the Q n distribution usually range 
from a full Q 4 (the silica network) to Q° network, the or- 
thosilicate glass, depending on the amount of soda x. At 
q = 2.5 g/cm 3 , the Q°, Q 5 and Q 6 populations were found 
to be negligible (less than 0.1% in each case), Q 1,2,3 ' 4 
populations from simulation being given in Table |l] and 
compared with results from a previous simulation^ 7 -, with 
results of NMR studies 41 and with results 17 from a ran- 
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FIG. 3: (Color online) Distribution of IV, V and Vl-fold coor- 
dinated silicon atoms (a) and of I, II and Ill-fold coordinated 
oxygen atoms (b) with respect to density. Sodium atoms are 
not taken into account in the enumeration of the neighbors, 
so that O 1 refer to NBOs. 



dom model proposed by Lacy 4 ^. First, we notice that our 
findings differ slightly with those obtained by Cormacfeii 
using the same potential. The origin may be due to 
the fact that the system has a different thermal history 
(the cooling methodology is slightly different although 
the cooling rate is the same). However, both simulations 
are consistent with the random model. Differences with 
experimental data are important and these shifts have 
been found even in very large-scale simulations^. They 
have been attributed to the fact that the high cooling rate 
used in simulations induces a structure with a higher ef- 
fective temperature 43 so that the simulated Q™ statistics 
is the one of a high temperature frozen liquid. 



TABLE I: Proportion of Q" populations at room density. 



Q" 


Present MD 


MD Cormack 17 


NMR 41 


Random model 4 ^ 


Q 1 


1.288 


1.857 


0.000 


2.985 


Q 2 


18.598 


15.571 


4.776 


16.716 


Q 3 


44.067 


49.000 


68.358 


41.493 


Q 4 


35.908 


33.571 


26.567 


37.910 



When density changes, we observe that the changes in 
Q™ populations are correlated with the change in O and 
Si coordination numbers (Fig. [3]), i.e. they take place 
only for g > 3 g/cm 3 . Q™ populations do not show any 
significant changes at low density. At g = 2.5 g/cm 3 , the 
Q 5 proportion starts to increase however and reaches a 
maximum at g = 4.0 g cm 3 . Again we notice a clear cor- 
relation between the Q 5 population and the proportion 
of Si v . The Q 6 proportion only starts to increase from 
g — 3.1 g /cm 3 , thus showing a behavior similar to the 
one of the proportion of Si VI . 



Bond-angle distributions 
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FIG. 4: (Color online) Distribution of Q" populations with 
increasing density. 



We now focus on the bond-angle distributions (BAD) 
and their variations with density, which have been 
shown to be extremely sensitive in other tetrahedral 



system; 



,38.44 



Even at ambient pressure, it allows to un- 



derstand how the basic structures of the glass connect to 
each other. The O-Si-0 BAD, which characterizes the Si 
tetrahedrons, is shown in Fig. [5^ for three selected den- 
sities (g = 2.5, 3.5 and 4.5 g/cm 3 ). As expected at the 
lowest density (g = 2.5 g/cm 3 ), the distribution is sharp 
and the average O-Si-0 angle is close to the ideal 109.5° 
tetrahedral angle (see also Fig. ^,). At intermediate den- 
sities however (g w 3.5 g/cm 3 ), the O-Si-0 angle displays 
now a bimodal distribution with a peak still located at 
109°, reminiscent of the initial tetrahedral structure, and 
a growing peak at 90° . The latter corresponds to the an- 
gle that is expected for an octahedral environment. An 
additional signature for this environment is provided by 
the contribution at 180° which appears for larger den- 
sities (Fig. [3^) and grows with g. At high density (g 
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FIG. 5: (Color online) O-Si-O (a), Si-BO-Si (b) and Si-NBO- 
Na (c) bond angle distributions at Q = 2.5, 3.5 and 4.5 g 
cm -3 . 



FIG. 6: (Color online) First (a) and second (b) moment of the 
O-Si-O, Si-BO-Si and Si-NBO-Na bond angle distributions. 



= 4.5 g/cm 3 , all silicon atoms display an octahedral en- 
vironment with a single peak at 90° (and the vanishing 
of the tetrahedral peak at 109°) and the contribution at 
180°. The second moment of the O-Si-0 BAD is shown 
in Fig. [SJa and grows from 5° at ordinary density (g = 2.5 
g/cm 3 ) up to more than 30° at high density, suggesting 
the appearance of a pressure induced disorder manifested 
by an increased angular excursion around a mean value. 
The Si-BO-Si angle characterizes the way two adjacent 
silicon tetrahedrons are connected. At ordinary density, 
the angle shows a broad distribution between 120° and 
180° and centered at 153° (see also Fig. [S^i), compared to 
the 142° experimental value from NMR 45 .^. This differ- 
ence was also observed by previous simulations, as re- 
viewed in^ and is consistent with the over estimated 
value of the Si-Si distance. Like the O-Si-0 angle, the 
Si-BO-Si angle displays a bimodal distribution as den- 
sity increases. The BAD shows a second peak close to 
100° at high density (g = 4.5 g/cm 3 ). This contribution 
is absent at intermediate densities (g = 3.5 g/cm 3 , see 
Fig. [5j}) but the trend with g is clearly correlated with 
the population of O m (see Fig. [5]d and Fig. [BJi). The 
decrease of the average value of the Si-BO-Si angle has 
also been observed experimentally in silicates 47 and for 



On the contrary, the Si-NBO-Na BAD does not show 
any significant change with density. 



6. Orientational parameter 

An interesting means to analyze the tetrahedral to oc- 
tahedral conversion in liquids and glasses is provided by 
the orientational order parameter q (introduced by Chau 
and Hardwick 4 - and rescaled in^i), which quantifies the 
extent to which a molecule and its four nearest neighbors 
adopt a tetrahedral arrangement. It is defined by : 



« = 1 -^E E (cos% fc + -)> (3) 

i=l k=j+l 

where 8%jk is the angle formed by the central Si atom 
i and its oxygen nearest neighbors j and k, the brack- 
ets representing an average over the central Si atoms i 
and over the time. This parameter is normalized so that 
its average varies between (randomly arranged bonds) 
and 1 (perfect tetrahedral network). It has been used 
for the analysis of diffusivity anomalies in relationship 
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FIG. 7: (Color online) q factor distributions for increasing 
densities. 
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FIG. 8: (Color online) First (left axis) and second (rigt axis) 
moment of the q factor distributions for increasing densities. 



with tetrahedral to octahedral changes in liquid water 34 , 
silica 33 and germania 4 ^. 

Fig. shows the calculated distribution of q values 
for 3 selected densities (g = 2.5, 3.5 and 4.5 g/cm 3 ). 
At room density, the distribution exhibits only a sharp 
peak close to q = 0.7 (see Fig. [8]), which corresponds 
to near-perfect tetrahedral order, as also found for silica 
at ambient pressure^ 3 -. At g = 3.5 g/cm 3 , the observed 
bimodal distribution suggests the existence of both tetra- 
hedral Si (high q peak) and higher coordinated Si (low q 
peak with q < 0.6). At g — 4.5 g/cm 3 , the tetrahedral 
contribution has vanished. The second moment of the 
distributions according to density is shown in Fig. |8] It 
characterizes the orientational disorder around central Si 
atoms, with tr q going from 0.02 at low density up to 0.1 



at g = 4.5 g/cm 3 , suggesting once again the appearance 
of a pressure induced disorder. 



B. Reciprocal space properties 

To investigate the structure of the glass on intermedi- 
ate length scales, the neutron structure factor has been 
computed. The partial structure factors have been first 
calculated from the pair distribution functions gij{r) : 



S l3 {Q) = 1 + eo 



fl 47rr 2 ( 3lJ (r)-l)^^F L (r)dr (4) 



where Q is the scattering vector, go is the average atom 
number density and R is the maximum value of the in- 
tegration in real space (here R — 15A). The Fi,(r) — 
sin(7rr / R) / (wr / R) term is a Lortch-type window function 
used to reduce the effect of the finite cutoff of r in the 
integration 50 . As discussed in^i, the use of this function 
reduces the ripples at low Q but induces a broadening of 
the structure factor peaks. The total neutron structure 
factor can then be evaluated from the partial structure 
factors following : 



Sn(Q) = (2J CiCjbibj) 1 2J CiCjhbjSijiQ) (5) 

where Cj is the fraction of i atoms (Si, O or Na) and bi 
is the neutron scattering length of the species (given by 
5.803, 4.1491 and 3.63 fm for oxygen, silicon and sodium 
atoms respectively 5 ^). 



1. Neutron structure factor 

The total neutron structure factor 5n for different in- 
creasing densities are shown on Fig. [§1 The room den- 
sity structure factor is compared both with the neutron 
diffraction results from Wright et ali^£ and the simulated 
NS2 glass from Horbach et alJ^, using an alternative 
(BKS) potential. We note that the agreement between 
simulation and experiment is good. The agreement of 
the first peak with experiment is discussed in details in 
the FSDP section below. The second peak position is well 
reproduced (3.0A _1 experimentally, compared to 3.0A _1 
from the present potential and 2.9A -1 from the BKS po- 
tential). The third peak position is also very well repro- 
duced (5.4A~ X experimentally, compared to 5.3A -1 from 
the present potential and 5.2A -1 from the BKS poten- 
tial). 

As observed on Fig. [SJ the density mainly influences 
the low wave vector part of the structure factors, which 
suggests that the main effects of density do not apply at 
short length scales. The shape in the high Q limit (at 
<5>10A _1 is nearly unchanged for g ^ 3.5 g/cm 3 . All 
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FIG. 9: (Color online) Neutron structure factors for increas- 
ing densities. Neutron diffraction results from Wright et 
al.— (open circles) and simulation results from Horbach at 
al.— (dotted line) are shown for comparison. Examples of 
Lorentzian fits of the FSDP are displayed in orange. 



peaks are shifted to higher wave vector (lower r) as den- 
sity increases, which is linked to the compaction of the 
network. Interestingly, the second moments of the differ- 
ent peaks do not show the same behavior with density. 
The so-called first sharp diffraction peak (FSDP) at very 
low Q becomes broader and less intense with increasing 
density, as discussed below. The main peak around 3 A -1 
becomes narrower as density increases, whereas the third 
one (around 5A _1 ) becomes broader. These trends can 
be analyzed in more details from the partial structure 
factors (see below). The shape of the other peaks are 
almost unaffected by density. 



2. Partial structure factors 

Fig. [TU1 shows the decomposition of the total structure 
factor into contributions of different pair structure factors 
Ssi-o(Q), So-o{Q) an d 5*Na-o( ( 5) f° r different increasing 
densities. The partial Ssi-Si(Q), <SNa-Na(Q) and Ssi-Na(Q) 
decay the fastest and have thus not been displayed. At 
normal density, the shape of these pair structure factors 
and the positions of the peaks are in excellent agreement 
with previously reported results from MD simulations 53 . 
At room pressure, the partial Ssi-o(Q) shows the most 
significant variations with Q both at short wave vector, 
correlated to the medium-range order of the silicate net- 
work, and at long Q, correlated to the strong short-range 
Si-0 order. As already observed in the total structure 
factor, most of the peaks are also shifted to higher Q 
(lower r) as density increases. 

The decomposition of the total structure factor can 
serve to understand the behavior of the main peak (~ 
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FIG. 10: (Color online) Partial structure factors Si-O, 0-0 
and Na-0 for increasing densities. 



3 A 1 ) and of the second main peak (~ 5 A 1 ) with den- 
sity. Indeed, the peak at 3A~ X in 0-0 and Na-0 par- 
tial structure factors becomes narrower as density in- 
creases, an effect which is related to the increased struc- 
tural medium-range order at high density. These peaks 
contribute the most of the second peaks of the total struc- 
ture factor. On the other hand, the main contribution 
for the peak at 5 A -1 of the total structure factor arises 
from the second Si-0 partial structure factor peak, which 
becomes broader as density increases. This apparent dis- 
order may be attributed to the appearance of coexisting 
tetrahedral and octahedral Si-0 environment as density 
increases. 



3. First sharp diffraction peak 

FSDPs are not simply the first of the many peaks 
of any diffraction pattern but display many anomalous 
behavior as a function of temperature, pressure and 
composition. 54 Since the position of the FSDP Qfsdp 
is smaller than Qp (the position of the principal peak of 
the structure factor, associated to the nearest-neighbor 
distance), the FSDP corresponds to structural correla- 
tions on a larger length scale. This feature has been ob- 
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FIG. 11: (Color online) FSDP position of the total Neutron 
structure factor and positions of each relevant partial struc- 
ture factors FSDP. The insert shows the associated character- 
istic distance d — 27t/Qfsdp of the total and partial structure 
factors. 



FIG. 12: (Color online) (a) Intensity of the FSDP of the total 
and partial structure factors, (b) FSDP FWHM of the total 
and partial structure factors. The insert shows the correlation 
length L = 2tt/FWHM. 



served both in covalen t 55 ' 56 and ionio^ 7 - amorphous sys- 
tem. In ionic systems, this medium range order has been 
associated to the forced separation between cations be- 
cause of their mutual Coulomb repulsion, thus producing 
a prepeak in the cation-cation structure factor—. Pre- 
peaks can also arise from size effects of the atoms of the 
network^. However, the network formation itself can 
have a major role since the FSDP is also observed in the 
monoatomic tetravalent systems a-Si and a-G o 60 i 61 . The 
FSDP origin is now usually explained by using a void- 
based model 5 -* 6 ^ in which ordering of interstitial voids 
occurs in the structure. 

The FSDPs we obtained from simulations were further 
studied by fitting them with Lorentzian functions (exam- 
ples of fitted functions can be seen on Fig. [9]). This choice 
is supported by the fact that the experimental results in 
neutron scattering factor of silica can be better fitted 
with a Lorentzian function than with a Gaussian one§2,. 
It should be noted that the fit has been done on the low 
Q part of the FSDP to avoid the contribution of the fol- 
lowing peaks. This allows to track precisely intensity, 
position and full-width at half maximum (FWHM) with 



density. Fig. [Til and [T^l show the position QfsdPj the 
intensity Jfsdp and the FWHM of the FSDP. The com- 
puted FSDP position at room density (l.SSA -1 ) is found 
to be in very good agreement with the one obtained from 
experiment (1.83A -1 )^ and with the one from the MD 
work by Corrales et al. (1.77A _1 )£i. Except at very low 
density (g < 2.5 g/cm 3 , negative pressure domain), the 
FSDP position increases with density while its intensity 
decreases. This trend is consistent with X-ray diffraction 
results from Benmore in densified silica. 6 — Interestingly, 
the FHWM of the FSDP exhibits a density window be- 
tween 2.3 and 3.3 g/cm 3 with a minimum found at g = 
2.7 g/cm 3 . 

Coming back to the real space correlations, the FSDP 
peak position Qfsdp is usually related to a character- 
istic repeat distance d = 27t/Qfsdp and the FWHM 
to a correlation length L = 27r/FWHM, sometimes 
also called 'coherence length', due to atomic density 
fluctuations^ 6 -^. The effect of irradiatio n 36 ! 68 , water 
conten t 67 ! 69 and alkali contend on the FSDP have been 
studied, leading to the idea that a depolymerization of 
the network (a decrease of the atomic order) is associ- 
ated to a decrease of the intensity of the FSDP and a 
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decrease in the characteristic distance d. A global un- 
derstanding of the correlation length L is lacking since it 
has been found to decrease with increasing potassium 
amount in silica network, to increase with increasing 
lithium amount and to show a maximum in sodium sili- 
cate when x — 0.2C&. It seems therefore highly system 
dependent. 

The inserts of Fig. [TT] and [T2b show the computed 
characteristic distance d and characteristic correlation 
length L as a function of density. It can be observed 
that the characteristic distance d decreases with density, 
which suggests a decrease of the medium range order 
(MRO). On the contrary, the correlation length L does 
not follow a general behavior with density since it shows 
a density window between 2.3 and 3.3 g/cm 3 with a max- 
imum at 2.7 g/cm 3 . We notice that the density of the 
maximum of L corresponds to the density of the begin- 
ning of the growth of the Si v fraction (see Fig. [jjt). 



4- Contributions to the FSDP 

Even though it can be noticed from Fig. QJj] that all 
partial structure factors show a FSDP, they do not con- 
tribute to the FSDP of the total structure factor at the 
same level. To understand the behavior of the FSDP, 
the position, intensity and FWHM of the FSDPs of each 
partial structure factor Sij(Q) have been computed. 

Fig. [TT] and [T2"h show the position and the inten- 
sity of the FSDPs according to density. At low den- 
sity (g < 2.7g/cm 3 ), the main contribution to the total 
FSDP clearly comes from the Si-0 FSDP, since their po- 
sition and intensity are similar and show the same trend. 
However, at larger densities, the partial FSDPs positions 
show a maximum (around g = 3.1g/cm 3 ) whereas the 
total FSDP continuously increases. These maximums 
correspond to minimums of the characteristic repeat dis- 
tance d and we notice that they occur at the density at 
which the Si VI fraction starts to grow. This shows that 
the total FSDP is not a simple superposition of the par- 
tial FSDP. The shift of the total FSDP to higher Q at 
high density can mainly be explained by the increased 
contribution of the main peak of the 5*o-o ( a t 3A _1 ) 
whose intensity grows with the density. 

Even though the link between the total and the par- 
tial FSDPs is not simple, it is interesting to notice that 
each partial FSDPs show a minimum of their FWHMs 
according to the density (see Fig. [T2"b). Si-0 and O- 
Na partial FSDPs FWHM reach their minimums around 
2.7g/cm 3 , corresponding to the minimum of the to- 
tal FSDP FWHM. The 0-0 partial FSDP shows the 
sharpest minimum of its FWHM around 3.2g/cm 3 (i.e. 
at larger density than for the total FSDP) but does not 
contribute a lot to the total FSDP due to its low intensity 
(see Fig. [TJi). 
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FIG. 13: (Color online) Vibrational density of states at room 
pressure computed from the Fourier transform of the velocity 
autocorrelation function (VAF) and from the diagonalization 
of the dynamical matrix (DM). The results are compared with 
the one from the simulation of Zotov et al. 70 using a different 
potential. The partial VDOS for Si, BO, NBO and Na are 
also shown. 



C. Vibrational properties 

The nature of the vibrational excitations of silicate 
glasses has so far remained a challenging issue. As con- 
trary to crystals, the lack of long-range structural or- 
der in amorphous solids strongly affects their vibrational 
dynamics. The appearance of an excess of vibrational 
modes over the Debye level at terahertz frequencies, the 
so-called Boson peak (BP), is one of the special features 
exhibited by glasses. 



1. Vibrational density of states 

The vibrational density of states (VDOS) g(oS) can be 
computed in two different ways. Starting from a relaxed 
glass (via energy minimization or cooling to OK), one 
can compute the dynamical matrix (DM) by evaluating 
the second derivative of the total energy with respect to 
small atomic displacements^. The diagonalization of the 
DM provides the eigenvalues, i.e. the frequency of each 
normal vibrational mode. Another way is to compute the 
Fourier transform of the velocity autocorrelation function 
(VAF) : 



I N fOO 

SM = Nh T ^ m j < v i W v j (0) > exp(icji) dt 

(6) 
where N is the number of atoms, rrij is the mass of an 
atom j, uj is the frequency and Vj(t) is the velocity of 
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FIG. 14: (Color online) Vibrational density of states, com- 
puted from the VAF, for different selected densities. 



FIG. 15: (Color online) Boson peak visualization for different 
selected densities. 



an atom j. It has been reported that both methods lead 
to quite similar VDOS in silica 72 . Although the DM ap- 
proach is far more expensive computationally than the 
VAF one, it should be noted that, looking at the eigen- 
vectors e^ associated to each eigenvalue frequency Wj, one 
can get the details of the nature of each normal mode. 
One can, for example, compute the partial VDOS g a (ui) 
(a = Si, BO, NBO, Na) for each atom defined as : 



t (ui) =g(ui) 



jGa 



/(w<)| S 



(7) 



where Bj(u>i) are the 3-component real space eigenvectors 
associated to the atoms a. 

Fig. [l~3l shows the VDOS, scaled to one, computed us- 
ing the two previously described methods. We observe a 
fair agreement between both methods, especially at low 
and intermediate frequency. The difference at high fre- 
quency can be explained by the harmonic assumption on 
which the DM method relies. The results are compared 
with the one from the simulation of Zotov et al^. We 
observe that the agreement is very poor even if some 
trends are similar (sharp peak at high frequency and ap- 
pearance of a new peak at low frequency increasing with 
respect to the amount of sodium) . We note that the sim- 
ulation from Zotov et al. uses a different potential (from 
Vessal) involving both 2- and 3-body terms and that the 
system is quite smaller (1080 atoms, as compared to 3000 
in the present simulation). Unfortunately, to our knowl- 
edge, no experimental VDOS is currently available for 
this composition. 

Using Eq. \7\ the partial VDOS have been computed 
and are shown on Fig. [13] Note that O atoms have 
been split into BOs and NBOs. Relying on the vi- 



brational analysis of silic; 



,72,73 



one can interpret some 



features of the present VDOS. Si atoms contribute the 



most at high frequency (27-37 THz, symmetric and anti- 
symmetric stretching modes) and at intermediate fre- 
quency (22 THz, O-Si-0 bending mode). BO atoms pre- 
dominant contributions occur at high frequency (29-37 
THz, symmetric stretching modes) and at low and inter- 
mediate frequency (0-26 THz, Si-BO-Si bending mode 
and symmetric stretching mode). The contribution of 
NBOs differs from the one of the BOs because of the 
decreased number of O-Si stretching modes and of the 
softening of the Si-NBO-Na bending mode as compared 
to the Si-BO-Si mode. Most of the low frequency contri- 
bution comes from Na atoms (0-10 THz, Na-NBO low- 
energy stretching modes). 

Fig. [141 shows the vibrational density of states, scaled 
to one, for different increasing densities. We observe 
some trends which are similar to the ones observed 
in densified silica 74 : decrease of the number of low- 
frequency modes, disappearance of the gap between in- 
termediate and high frequency modes around 27 THz 
and broadening of the high-frequency peak. The low- 
frequency region, whose contribution mainly comes from 
Na atoms, is the most affected part of the VDOS, suggest- 
ing that the Na vibrational modes are strongly modified 
during densification. 

As density increases, the low-frequency peak coming 
from Na-0 bounds decreases in intensity and is shifted 
to higher frequencies. The high frequency peak coming 
from Si-0 bounds becomes broader but does not show 
any significant frequency shift. 



2. Boson peak 

The origin of the BP in silica, still controversial, has 
been associated to the existence of local modes involving 
rocking motions of distorted Si04 tetrahedrons 7 ^— . The 
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FIG. 16: (Color online) Boson peak position (left axis) and 
intensity (right axis) with respect to density. 



BP can be observed by looking at the excess VDOS over 
the Debye law g(u>) — 9n(u) or at the quantity g(uj)/u> . 
The latter quantity can be identified with the one-phonon 
scattering cross section as measured in neutron scatter- 
ing experiment a 77 ' 78 and is shown on Fig. [15] for three 
selected densities. A pronounced peak can be observed 
at each density, even if its intensity decreases with den- 
sity. At room pressure, the BP is found to be located 
at wbp=1-3 THz, which is lower than the value found 
experimentally (Raman) of 1.95 THz (65 cm _1 )2. 

The BP properties have been further analyzed by com- 
puting its position wbp and its intensity Ibp, quantities 
that are displayed on Fig. [TBI It can been observed that 
Ibp decreases with density, while wbp increases with den- 
sity. Both of these two trends (decrease of the intensity 
and increase of the frequency) have been observed ex- 
perimentally in many system, such as in pure silica^, in 
lithium silicate glass^i, in a Na2FeSi30s glass^ as well 
as in different polymers 83 . 



IV. LIQUID 



A. Thermodynamics 



FIG. 17: (Color online) Isotherms for glass and liquid NS2. 
The curves are separated by 500K each. The inset shows the 
corresponding data in (g, P) together with the BM fits (solid 
lines) . 
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FIG. 18: (Color online) Isothermal compressibility kt with 
respect to density in liquid NS2 for various temperatures rang- 
ing from 3000 to 1500K separated by 500K each. The curves 
are computed using the BM EOS. The inset shows the room 
pressure density go change with temperature. 



To evaluate the equation of state (EOS) of the liquid, 
many thermodynamics points have been computed (T, g, 
P). The following range has been studied : 1.5 ^ g ^ 5.5 
g/cm 3 and 1500 ^J T ^J 3000 K, which correspond to the 
following pressure range : —2.23 ^ P ^ 150 GPa. In con- 
trast with previous works on molecular fluids 84,85 and sil- 
ica, where the data were fitted using a Van der Waals type 
EOS, the data of the current simulations were fitted with 
a Birch-Murnaghan equation of state (BM EOS) that has 
a simpler for m 86 ' 87 . It has revealed to give reasonable fits 
in the case of a liquid densified germania— and is widely 



used in geophysical studies (see for example 89 ). 

Fig. IT7l shows the isotherms of the glass and the liquid 
in the (P, V) representation. The data have been fitted 
far from the critical region with the BM EOS, that has 
the following form : 



P = |^((^) 7 / 3 - [if /3) (1 .J (4 _ Kl ){[^ 1)) 
2 go go 4 g 

(8) 
where K is the bulk modulus at P=0, K\ — dK/dP at P 
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FIG. 19: (Color online) Radial distribution function for in- 
creasing temperatures. 



= and go is zero-pressure density of the liquid. The fit 
can be made with two parameters only (K and K{) since 
Qo can be accessed from the isothermal data displayed on 
Fig. [lTj It can be seen that the data are very well fitted 
by the BM EOS along all the density and temperature 
range. 

In addition, the BM EOS allows to have access to the 
bulk modulus K at P=0 and to the isothermal compress- 
ibility kt = g~ 1 (dg/dP)T according to the density (plot- 
ted on Fig. [T8]) . The observed behavior, enhanced com- 
pressibility with falling density, is realistic, as well as the 
decrease of the bulk modulus at P = with respect to 
the temperature (see the insert of Fig. \TE\i . The results 
are in good agreement with the only experimental data 
on liquid NS2 we are aware of (K = 13.4 GPa and kt = 
0.075 GPa" 1 at T=1500 K and P = Q~p.. 



FIG. 20: (Color online) Distribution of IV, V and Vl-fold 
coordinated silicon atoms (a) and of I, II and Ill-fold coordi- 
nated oxygen atoms (b) with respect to density both at 300K 
(filled symbols) and 2000K (open symbols). Sodium atoms 
are not taken into account in the enumeration of the neigh- 
bors, so that O 1 refer to NBOs. 



icant shift as temperature increases. The only behavior 
that can be seen is a broadening of all peaks as tempera- 
ture increases, which can be explained by the increasing 
disorder due to the increasing thermal energy. The same 
trends are observed in pure silica^. 

Eventually, the influence of the temperature on the co- 
ordination numbers has been checked. The populations 
of the different Si and O species according to the density 
are plotted on Fig. [Ulboth at T = 300K and 2000K. The 
transition between Si IV into Si v and Si VI is still clearly 
observed, although the transition occurs at a lower den- 
sity than in the glass (density shift of approximately 0.2 
g/cm 3 ). The same shift can be observed for the O species. 



B. Structure 

Fig. [ll?l shows the total correlation function g^{r) for 
different increasing temperatures, both at low and high 
density. It can be observed that density seems to have 
a more critical influence on structure than temperature. 
Indeed, the positions of the peaks do not show any signif- 



V. CONCLUSION 

Our purpose in the present paper has been to provide a 
systematic and extensive study of the properties of den- 
sified glassy and liquid NS2 sodium silicate. 

While bond distances remains nearly unchanged, pres- 
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sure has a strong effect on angles and coordination num- 
bers. A transition from tetrahedral to octahedral silicon 
environment is found. The fraction of NBOs decreases 
and O m tricluster are observed at high density. The 
usual vibrational behavior is observed, i.e. the decrease 
of the amount of low frequency modes, the increase of the 
frequency of the Boson peak and the decrease of its in- 
tensity under pressure. Expected anomalous effects are 
found in the medium range order (increase of the po- 
sition of the FSDP and decrease of its intensity under 
pressure), but, more surprisingly, we observe a minimum 
of the FWHM of the FSDP according to the density. 
Temperature is found to have only small effects on struc- 
ture and the Birch-Murnaghan equation of state allows 
to reproduce the densification of the liquid at each tem- 
perature. 

Finally, it is worth mentioning that, as ambient pres- 



sure, it is well-known that changes in composition of the 
glass (and especially in the amount of sodium atoms) in- 
duce changes of the degree of polymerization of the glass. 
These competitive effects (depolymerization by sodium 
atoms and repolymerization by the pressure) should be 
addressed in the future for a better understanding of the 
glass network properties. 
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